#include <UTILS/SBuffer.h>
#include <PECAgent/PEContext.h>
#include <DPS/PH/DpsFunctions.h>
#include <DPS/PH/DpsOPContext.h>
#include "UDOP_OpGlobalX.h"
#include "analysis.h"

SPC_NAMESPACE_USE;
UTILS_NAMESPACE_USE;
using namespace DPS;
using namespace std;
using namespace estd;

/*
 * Function used to initialize the operator
 */
void UDOP_OpGlobalX::initializeOperator() 
{
  // Get the paramter's value
  char str10[150];
  const DpsOPContext &context = this->getContext();
  const std::string &args = context.getCmdArguments();
  strcpy(str10, args.c_str());
  sscanf(str10,"%*s %*c%d%*c %*s %*c%d%*c",&num_ensemble,&astep);
  
  // Allocate memory space
  //double **tmpmass;
  tmpmass = (double **)malloc(num_ensemble*sizeof(double *));
  //tmpmass = new double*[num_ensemble];
  
  //float **tmptlat;
  tmptlat = (float **)malloc(num_ensemble*sizeof(float *));
  //tmptlat = new float*[num_ensemble];
  
  //float **tmptlon;
  tmptlon = (float **)malloc(num_ensemble*sizeof(float *));
  //tmptlon = new float*[num_ensemble];
  
  //float **tmpzlvl;
  tmpzlvl = (float **)malloc(num_ensemble*sizeof(float *));
  //tmpzlvl = new float*[num_ensemble];

  //float **tmpsigh;
  tmpsigh = (float **)malloc(num_ensemble*sizeof(float *));
  //tmpsigh = new float*[num_ensemble];
  
  //float **tmpsigw;
  tmpsigw = (float **)malloc(num_ensemble*sizeof(float *));
  //tmpsigv = new float*[num_ensemble];

  //float **tmpsigv;
  tmpsigv = (float **)malloc(num_ensemble*sizeof(float *));
  //tmpsigv = new float*[num_ensemble];

  //int16_t **tmppage;
  tmppage = (int16_t **)malloc(num_ensemble*sizeof(int16_t *));
  //tmppage = new int16_t*[num_ensemble];

  //int16_t **tmphdwp;
  tmphdwp = (int16_t **)malloc(num_ensemble*sizeof(int16_t *));
  //tmphdwp = new int16_t*[num_ensemble];

  //int16_t **tmpptyp;
  tmpptyp = (int16_t **)malloc(num_ensemble*sizeof(int16_t *));
  //tmpptyp = new int16_t*[num_ensemble];

  //int16_t **tmppgrd;
  tmppgrd = (int16_t **)malloc(num_ensemble*sizeof(int16_t *));
  //tmppgrd = new int16_t*[num_ensemble];

  //int *en_numpar;
  en_numpar = (int *)malloc(num_ensemble*sizeof(int));
  
  int i = 0;
  for(i=0; i<num_ensemble; i++)
  {
	  tmpmass[i] = NULL;
	  tmptlat[i] = NULL;
	  tmptlon[i] = NULL;
	  tmpzlvl[i] = NULL;
	  tmpsigh[i] = NULL;
	  tmpsigw[i] = NULL;
	  tmpsigv[i] = NULL;
	  tmppage[i] = NULL;
	  tmphdwp[i] = NULL;
	  tmpptyp[i] = NULL;
	  tmppgrd[i] = NULL;
  }
  //en_numpar = new int[num_ensemble];
  // Create a user defined metric
  // opm.createMetric("intuples","Number of ingested tuples");
}

/*
 * Function used to finalize the operator
 */
void UDOP_OpGlobalX::finalizeOperator() 
{
	// free space
	int i = 0;
	for (i=0; i<num_ensemble; i++)
		{
		if (tmpmass[i] != NULL)
			free(tmpmass[i]);
			//delete[] tmpmass[i];
		
		if (tmptlat[i] != NULL)
			free(tmptlat[i]);
			//delete[] tmptlat[i];
		
		if (tmptlon[i] != NULL)
			free(tmptlon[i]);
			//delete[] tmptlon[i];
			
		if (tmpzlvl[i] != NULL)
			free(tmpzlvl[i]);
			//delete[] tmpzlvl[i];
		
		if (tmpsigh[i] != NULL)
			free(tmpsigh[i]);
			//delete[] tmpsigh[i];
		
		if (tmpsigw[i] != NULL)
			free(tmpsigw[i]);
			//delete[] tmpsigw[i];
		
		if (tmpsigv[i] != NULL)
			free(tmpsigv[i]);
			//delete[] tmpsigv[i];
		
		if (tmppage[i] != NULL)
			free(tmppage[i]);
			//delete[] tmppage[i];
		
		if (tmphdwp[i] != NULL) 
			free(tmphdwp[i]);
			//delete[] tmphdwp[i];
		
		if (tmpptyp[i] != NULL)	
			free(tmpptyp[i]);
			//delete[] tmpptyp[i];
		
		if (tmppgrd[i] != NULL)
			free(tmppgrd[i]);
			//delete[] tmppgrd[i];
		}

	free(tmpmass); free(tmptlat);
	free(tmptlon); free(tmpzlvl);
	free(tmpsigh); free(tmpsigw);
	free(tmpsigv); free(tmppage);
	free(tmphdwp); free(tmpptyp);
	free(tmppgrd); free(en_numpar);

}

/*
 * Function used to process command line arguments
 */
void UDOP_OpGlobalX::processCmdArgs(const string & /*args*/) 
{
  SPCDBG(L_INFO, "Process command line arguments...", "dpsop");
  // ArgContainer arg(args); // include "ArgContainer.h"
  // arg.argc and arg.argv are now available for use
  // do other initialization here
}

/*
 * Process function for input port 0
 */
void UDOP_OpGlobalX::processInput0(const IPort0_t & tuple)
{
  // declare for output tuples for the averaged ones
  const dlist<double> & kmass=(const dlist<double>&) (tuple.get_mass());
  const dlist<float> & ktlat =(const dlist<float>&) (tuple.get_tlat());
  const dlist<float> & ktlon =(const dlist<float>&) (tuple.get_tlon());
  const dlist<float> & kzlvl =(const dlist<float>&) (tuple.get_zlvl());
  const dlist<float> & ksigh =(const dlist<float>&) (tuple.get_sigh());
  const dlist<float> & ksigw =(const dlist<float>&) (tuple.get_sigw());
  const dlist<float> & ksigv =(const dlist<float>&) (tuple.get_sigv());
  const dlist<int16_t> & kpage =(const dlist<int16_t>&) (tuple.get_page());
  const dlist<int16_t> & khdwp =(const dlist<int16_t>&) (tuple.get_hdwp());
  const dlist<int16_t> & kptyp =(const dlist<int16_t>&) (tuple.get_ptyp());
  const dlist<int16_t> & kpgrd =(const dlist<int16_t>&) (tuple.get_pgrd());
  int kstreamid = tuple.get_stid();
  
  // Obtain the values of input tuples 
  
  int i = 0;
  int streamid = 0;
  
  // Obtain stream from all input ports
  //otuple->assign(*ituple);
  if (kmass.size()>0)
	{
    streamid = kstreamid;
    en_numpar[streamid] = kmass.size();
    
  tmpmass[streamid]=(double *)realloc(tmpmass[streamid],
  					en_numpar[streamid]*sizeof(double));
  					
  tmptlat[streamid]=(float *)realloc(tmptlat[streamid],
  					en_numpar[streamid]*sizeof(float));
  										
  tmptlon[streamid]=(float *)realloc(tmptlon[streamid],
  					en_numpar[streamid]*sizeof(float));
  
  tmpzlvl[streamid]=(float *)realloc(tmpzlvl[streamid],
  					en_numpar[streamid]*sizeof(float));
  
  tmpsigh[streamid]=(float *)realloc(tmpsigh[streamid],
  					en_numpar[streamid]*sizeof(float));
  										
  tmpsigv[streamid]=(float *)realloc(tmpsigv[streamid],
  					en_numpar[streamid]*sizeof(float));
  					
  tmpsigw[streamid]=(float *)realloc(tmpsigw[streamid],
  					en_numpar[streamid]*sizeof(float));
  					
  tmppage[streamid]=(int16_t *)realloc(tmppage[streamid],
  					en_numpar[streamid]*sizeof(int16_t));
  					
  tmphdwp[streamid]=(int16_t *)realloc(tmphdwp[streamid],
  					en_numpar[streamid]*sizeof(int16_t));
  					
  tmpptyp[streamid]=(int16_t *)realloc(tmpptyp[streamid],
  					en_numpar[streamid]*sizeof(int16_t));
  
  tmppgrd[streamid]=(int16_t *)realloc(tmppgrd[streamid],
  					en_numpar[streamid]*sizeof(int16_t));
  					
  for(i=0;i<en_numpar[streamid];i++)
     {
     tmpmass[streamid][i]=kmass[i];
	 tmptlat[streamid][i]=ktlat[i];
	 tmptlon[streamid][i]=ktlon[i];
	 tmpzlvl[streamid][i]=kzlvl[i];   
     tmpsigh[streamid][i]=ksigh[i];
     tmpsigw[streamid][i]=ksigw[i];
	 tmptlon[streamid][i]=ksigv[i];
	 tmppage[streamid][i]=kpage[i];
	 tmphdwp[streamid][i]=khdwp[i];
	 tmpptyp[streamid][i]=kptyp[i];
	 tmppgrd[streamid][i]=kpgrd[i];
	 }

  gflag = true;	 
	 
	} // end of if (array size not zero) 
  else
	{
		gflag = false; 
    	streamid = -1;
	}

	cnt++;

	if(cnt%num_ensemble == 0 )
	{
		handleWindow();
	}

}

/*-------------------------------------------------------*/
/*-------------------------------------------------------*/
void UDOP_OpGlobalX::handleWindow()
{
  //1.find the 1st longest numpar
  int i = 0;
  int temp=0; int kernel = 0;
  for (i=0; i<num_ensemble; i++)
  {
	if ( en_numpar[i]> temp)
	{
		temp = en_numpar[i];
		kernel = i;
	}
  }
 
  
//  switch (value)
//  {
//  case Punctuation::WINDOW_MARKER:
  		if (gflag)// && (tstep%astep==0))
		{
		BackgroundEn(temp, kernel);
//		BackgroundEn(temp, kernel);			
		}
//  		break;
//  case Punctuation::NON_EXISTENT:
//  		break;
//  case Punctuation::INVALID:
//  		break;					
//  }


	
}

void UDOP_OpGlobalX::BackgroundEn(int temp, int kernel)
//void UDOP_OpGlobal::BackgroundEn(int temp, int kernel)		
{
  // assign the kernel-th member's value to it
  int i, j;
 //2.Based on this longest PARDUMP, allocate space
  dlist<double> xb_bar; xb_bar.reserve(temp);
  dlist<float> xb_tlat; xb_tlat.reserve(temp);
  dlist<float> xb_tlon; xb_tlon.reserve(temp);
  dlist<float> xb_zlvl; xb_zlvl.reserve(temp);
  dlist<float> xb_sigh; xb_sigh.reserve(temp);
  dlist<float> xb_sigw; xb_sigw.reserve(temp);
  dlist<float> xb_sigv; xb_sigv.reserve(temp);
  dlist<int16_t> xb_page; xb_page.reserve(temp);
  dlist<int16_t> xb_hdwp; xb_hdwp.reserve(temp);
  dlist<int16_t> xb_ptyp; xb_ptyp.reserve(temp);
  dlist<int16_t> xb_pgrd; xb_pgrd.reserve(temp);
  dlist<int16_t> xb_idx; xb_idx.reserve(temp);
  
  //3. Allocate space for the rest
  dmatrix<double> en_mass = DPS::makeDoubleMatrix(num_ensemble,temp);
  dmatrix<float> en_tlat = DPS::makeFloatMatrix(num_ensemble,temp);
  dmatrix<float> en_tlon = DPS::makeFloatMatrix(num_ensemble,temp);
  dmatrix<float> en_zlvl = DPS::makeFloatMatrix(num_ensemble,temp);
  dmatrix<float> en_sigh = DPS::makeFloatMatrix(num_ensemble,temp);
  dmatrix<float> en_sigw = DPS::makeFloatMatrix(num_ensemble,temp);
  dmatrix<float> en_sigv = DPS::makeFloatMatrix(num_ensemble,temp);
  dmatrix<int16_t> en_page = DPS::makeShortMatrix(num_ensemble,temp);
  dmatrix<int16_t> en_hdwp = DPS::makeShortMatrix(num_ensemble,temp);
  dmatrix<int16_t> en_ptyp = DPS::makeShortMatrix(num_ensemble,temp);
  dmatrix<int16_t> en_pgrd = DPS::makeShortMatrix(num_ensemble,temp);

  for (i=0; i<temp; i++)
  {
  xb_bar[i] = tmpmass[kernel][i];
  en_mass(kernel,i) = xb_bar[i];
  xb_tlat[i]= tmptlat[kernel][i];
  xb_tlon[i]= tmptlon[kernel][i];
  xb_zlvl[i]= tmpzlvl[kernel][i];
  xb_sigh[i]= tmpsigh[kernel][i];
  xb_sigw[i]= tmpsigw[kernel][i];
  xb_sigv[i]= tmpsigv[kernel][i];
  xb_page[i]= tmppage[kernel][i];
  xb_hdwp[i]= tmphdwp[kernel][i];
  xb_ptyp[i]= tmpptyp[kernel][i];
  xb_pgrd[i]= tmppgrd[kernel][i];
  xb_idx[i] = 1;
  }
  
  //4. merge the rest members
  float tmpdist = 0.0;
  float dist =0.0;
  int idx = 0;
  int k = 0;
	for (i=0; i<num_ensemble; i++)
	{
	if(i!=kernel)
	{
		for (j=0; j<en_numpar[i]; j++)
		{	tmpdist = 9999999999.0; //reset;
			for (k=0; k<temp; k++)
			{
			dist =(tmptlat[i][j]-tmptlat[kernel][k]) * 
				(tmptlat[i][j]-tmptlat[kernel][k]) +
				(tmptlon[i][j]-tmptlon[kernel][k]) *
				(tmptlon[i][j]-tmptlon[kernel][k]) +
				(tmpzlvl[i][j]-tmpzlvl[kernel][k]) *
				(tmpzlvl[i][j]-tmpzlvl[kernel][k]);
			if (dist < tmpdist)
				{
				tmpdist = dist;
				idx = k;
				}
			} // end of k
			
	//append j-th of rest member to the idx-th element in kernel
			xb_bar[idx] += tmpmass[i][j];
			en_mass(i,idx) += tmpmass[i][j];
			xb_tlat[idx] += tmptlat[i][j];
			xb_tlon[idx] += tmptlon[i][j];
			xb_zlvl[idx] += tmpzlvl[i][j];
			xb_sigh[idx] += tmpsigh[i][j];
			xb_sigw[idx] += tmpsigw[i][j];
			xb_sigv[idx] += tmpsigv[i][j];
			xb_page[idx] += tmppage[i][j];
			xb_idx[idx] += 1;
			en_tlat(i,idx) = tmptlat[i][j];
			en_tlon(i,idx) = tmptlon[i][j];
			en_zlvl(i,idx) = tmpzlvl[i][j];
			en_sigh(i,idx) = tmpsigh[i][j];
			en_sigw(i,idx) = tmpsigw[i][j];
			en_sigv(i,idx) = tmpsigv[i][j];
			en_page(i,idx) = tmppage[i][j];
			en_hdwp(i,idx) = tmphdwp[i][j];
			en_ptyp(i,idx) = tmpptyp[i][j];
			en_pgrd(i,idx) = tmppgrd[i][j];
		} //end of j
	} // end of i
	else
	{
		//do nothing
	}
	}// end of num_ensemble

    //5. Average the result from step4. Notice. mass is divided by #ensemble
	//   the rests are divided by its according counter.
	for (i=0; i<temp; i++)
	{
		xb_bar[i] /= (double)xb_idx[i];
		xb_tlat[i] /= (float)xb_idx[i];
		xb_tlon[i] /= (float)xb_idx[i];
		xb_zlvl[i] /= (float)xb_idx[i];
		xb_sigh[i] /= (float)xb_idx[i];
		xb_sigw[i] /= (float)xb_idx[i];
		xb_sigv[i] /= (float)xb_idx[i];
		xb_page[i] /= (int16_t)xb_idx[i];
	}
	// Also the rest of en_pardump
	for (i=0; i<num_ensemble; i++)
	{
	if (i!=kernel)
	{
		for (j=0; j<temp; j++)
		{
		if (en_page(i,j) == 0)
			{
			en_mass(i,j) = xb_bar[j];
			en_tlat(i,j) = xb_tlat[j];
			en_tlon(i,j) = xb_tlon[j];
			en_zlvl(i,j) = xb_zlvl[j];
			en_sigh(i,j) = xb_sigh[j];
			en_sigw(i,j) = xb_sigw[j];
			en_sigv(i,j) = xb_sigv[j];
			en_page(i,j) = xb_page[j];
			en_hdwp(i,j) = xb_hdwp[j];
			en_ptyp(i,j) = xb_ptyp[j];
			en_pgrd(i,j) = xb_pgrd[j];
			}
			
		}
	}
	}
	//6. Get Xb_g = H_en - {xb_bar}
	for (i=0; i<num_ensemble; i++)
		for (j=0; j<temp; j++)
		{
		en_mass(i,j) -= xb_bar[j];
		}
		
		tstep++;
	/*	double test0 = en_mass(0,2);
		double test1 = en_mass(1,2);
		double test2 = en_mass(2,2);
		double test3 = en_mass(3,2);
		double test4 = en_mass(4,2);
		double test5 = en_mass(5,2);
		double test6 = en_mass(6,2);
		printf("tetst %lf %lf %lf %lf %lf %lf %lf\n", test0, test1, test2, test3, test4, test5, test6);*/
	OPort0_t otuple;
	otuple.set_xb_bar(xb_bar);
	otuple.set_xb_tlat(xb_tlat);
	otuple.set_xb_tlon(xb_tlon);
	otuple.set_xb_zlvl(xb_zlvl);
	otuple.set_xb_sigh(xb_sigh);
	otuple.set_xb_sigw(xb_sigw);
	otuple.set_xb_sigv(xb_sigv);
	otuple.set_xb_page(xb_page);
	otuple.set_xb_hdwp(xb_hdwp);
	otuple.set_xb_ptyp(xb_ptyp);
	otuple.set_xb_pgrd(xb_pgrd);
	otuple.set_en_mass(en_mass);
	otuple.set_en_tlat(en_tlat);
	otuple.set_en_tlon(en_tlon);
	otuple.set_en_zlvl(en_zlvl);
	otuple.set_en_sigh(en_sigh);
	otuple.set_en_sigw(en_sigw);
	otuple.set_en_sigv(en_sigv);
	otuple.set_en_page(en_page);
	otuple.set_en_hdwp(en_hdwp);
	otuple.set_en_ptyp(en_ptyp);
	otuple.set_en_pgrd(en_pgrd);
	otuple.set_flag(gflag);
	submit0(otuple);	
}
